Load obs and covar data

#Load Observed NYC case data
Observed_data = read.csv(paste0(
  "../Generated_Data/observed_data_",
  model_name, ".csv"))
head(Observed_data)

##   Y times obs_prop_positive
## 1 0     1                NA
## 2 0     2        0.00000000
## 3 2     3        0.25000000
## 4 2     4        0.05555556
## 5 7     5        0.15555556
## 6 0     6        0.00000000

### Define start date
true_start_date = as.Date("2020-03-01")
t0 = 0
start_of_year = as.Date("2020-01-01")
first_saturday_in_year = as.Date("2020-01-04")

## Compartment/Queue Cohort Numbers
M = 5
V = 13
K = 14


#Declare Csnippets and data
source("Csnippet_nyc_coronavirus_model_N_12.R")

## Load NYC covariate data
covariate_df = read.csv(file =
                          paste0("../Generated_Data/covariate_data_",
                                 model_name, ".csv"))



### Create covariate table
covar=covariate_table(
  time=covariate_df$times,
  L_advanced_2_days=covariate_df$L_advanced_2_days,
  F_w_y = covariate_df$F_w_y,
  L_orig = covariate_df$L_orig,
  w = covariate_df$Week,
  y = covariate_df$Year,
  times="time"
)

Plot 6 individual trajectories from big b_a parameter combination

select_trajectories  = filter(sim_data_big_b_a, .id %in% seq(from = 5, to = 10))
select_trajectories = dplyr::select(select_trajectories, time, .id, Y)
select_trajectories$type = "Sim"

library(RColorBrewer)
full_blue_pallete = brewer.pal(9, "Blues")
sim_traj_pallete = full_blue_pallete[9:4]
sim_traj_scale = scale_color_manual(name = "Legend", values = c( sim_traj_pallete), labels =  c("Sim_Traj_1", "Sim_Traj_2", "Sim_Traj_3", "Sim_Traj_4", "Sim_Traj_5", "Sim_Traj_6")) 
p = ggplot(data = select_trajectories,
           aes(x = time, y = Y, color = .id)) + geom_point(size = 2) + geom_line(aes(group = .id)) + rahul_theme + rahul_man_figure_theme + theme_white_background +
  sim_traj_scale
p

png("../Figures/Representative_Simulations/big_b_a_5_traj_sim.png")
print(p)
dev.off()

## quartz_off_screen 
##                 2

p = ggplot(data = select_trajectories,
           aes(x = time, y = Y, color = .id)) + geom_point(size = 2) + geom_line(aes(group = .id)) + rahul_theme +
  rahul_man_figure_theme + theme_white_background +
  sim_traj_scale + geom_point(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  geom_line(data = Observed_data,
                              aes(x = times, y = Y), color = 'red')
p

png("../Figures/Representative_Simulations/big_b_a_5_traj_sim_with_obs.png")
print(p)
dev.off()

## quartz_off_screen 
##                 2

p = ggplot(data = select_trajectories,
           aes(x = time, y = Y, color = .id)) + geom_point(size = 2) + geom_line(aes(group = .id))  +
  theme_white_background +
  sim_traj_scale +
    facet_wrap(~.id, ncol = 1) +
  geom_point(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  geom_line(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  theme(legend.position = "None") +
  theme(strip.background = element_blank(),
        strip.text.x = element_blank())
p

png("../Figures/Representative_Simulations/big_b_a_5_traj_sim_with_obs_facet.png")
print(p)
dev.off()

## quartz_off_screen 
##                 2

p = ggplot(data = select_trajectories,
           aes(x = time, y = Y, color = .id)) + geom_point(size = 2) + geom_line(aes(group = .id))  +
  theme_white_background +
  sim_traj_scale +
    facet_wrap(~.id, ncol = 1) +
  geom_point(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  geom_line(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  theme(legend.position = "None") 
p

Pick “good” simluation

(One that corresponds well with the data by eye)

big_b_a_single_traj = select_trajectories %>%
  filter(.id == "10")
p = ggplot(data = big_b_a_single_traj,
           aes(x = time, y = Y)) + geom_point(size = 2,
                                              color = 'blue') +
  geom_line(color = 'blue') + rahul_man_figure_theme +
  geom_point(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  geom_line(data = Observed_data,
                              aes(x = times, y = Y), color = 'red')
p

png("../Figures/Representative_Simulations/big_b_a_best_traj_vs_obs.png")
print(p)
dev.off()

## quartz_off_screen 
##                 2

Save single trajectory for high b_a parameter combination

NYC_full_testing_data =
  read.csv("../Generated_Data/NYC_full_testing_data.csv")
head(NYC_full_testing_data)

##    Test.Date New_Positives Cumulative_Number_of_Positives
## 1 03/02/2020             0                              0
## 2 03/03/2020             0                              0
## 3 03/04/2020             2                              2
## 4 03/05/2020             2                              4
## 5 03/06/2020             7                             11
## 6 03/07/2020             0                             11
##   Total_Number_of_Tests_Performed Cumulative_Number_of_Tests_Performed
## 1                               0                                    0
## 2                               8                                    8
## 3                               8                                   16
## 4                              36                                   52
## 5                              45                                   97
## 6                              64                                  161
##     Prop_Pos Not_Pos       Date
## 1         NA       0 2020-03-02
## 2 0.00000000       8 2020-03-03
## 3 0.25000000       6 2020-03-04
## 4 0.05555556      34 2020-03-05
## 5 0.15555556      38 2020-03-06
## 6 0.00000000      64 2020-03-07

nyc_rel_testing_data = NYC_full_testing_data %>%
  mutate(times = as.numeric(as.Date(Date) - true_start_date)) %>%
  select(times,
         Total_Number_of_Tests_Performed =
           Total_Number_of_Tests_Performed)

big_b_a_single_traj_data = big_b_a_single_traj %>%
  dplyr::select(Y =Y, times = time) %>%
  join(nyc_rel_testing_data) %>%
  mutate(obs_prop_positive = Y/Total_Number_of_Tests_Performed) %>%
  dplyr::select(-Total_Number_of_Tests_Performed)

## Joining by: times

big_b_a_single_traj_data$obs_prop_positive[1] = NA
big_b_a_single_traj_data

##       Y times obs_prop_positive
## 1     0     1                NA
## 2     0     2        0.00000000
## 3     1     3        0.12500000
## 4     0     4        0.00000000
## 5     3     5        0.06666667
## 6     8     6        0.12500000
## 7    11     7        0.17460317
## 8    15     8        0.11538462
## 9    30     9        0.16949153
## 10   22    10        0.12643678
## 11   26    11        0.07008086
## 12   54    12        0.07814761
## 13  182    13        0.35616438
## 14  132    14        0.16019417
## 15  289    15        0.32435466
## 16  310    16        0.14065336
## 17 1566    17        0.46830144
## 18 1300    18        0.26837325
## 19 1195    19        0.23537522
## 20 1363    20        0.19676628
## 21 3179    21        0.48064711
## 22 3733    22        0.69296454
## 23 5408    23        0.95110798
## 24 4619    24        0.66412653
## 25 3283    25        0.47421638
## 26 3243    26        0.40802718
## 27 3211    27        0.44584838
## 28 6236    28        1.00564425
## 29 6012    29        0.60776385
## 30 4645    30        0.63361069
## 31 5537    31        0.73698922
## 32 6506    32        0.69486276
## 33 7563    33        0.67238620
## 34 4923    34        0.63002304
## 35 6839    35        0.84089512
## 36 6747    36        0.75309744
## 37 4989    37        0.47190692
## 38 5744    38        0.53702319
## 39 2732    39        0.24775551
## 40 4991    40        0.42433260
## 41 4069    41        0.37735324
## 42 2762    42        0.34297777
## 43 4502    43        0.51687715
## 44 5301    44        0.32702036
## 45 4605    45        0.36769403
## 46 3848    46        0.32732222
## 47 4428    47        0.38005321
## 48 3402    48        0.33122383
## 49 5659    49        0.72569890
## 50 3178    50        0.35995016
## 51 2813    51        0.24993336
## 52 3722    52        0.27850943
## 53 2576    53        0.14220259
## 54 2383    54        0.13824110
## 55 2321    55        0.17111472
## 56 2826    56        0.26809601
## 57 2370    57        0.23958755
## 58 2811    58        0.19982939
## 59 1143    59        0.08014304
## 60 1931    60        0.14310064
## 61 1485    61        0.09026258
## 62 1254    62        0.09386930
## 63 1085    63        0.11815311
## 64 2151    64        0.19885366
## 65 1102    65        0.08447034
## 66 2190    66        0.13378948
## 67 1700    67        0.11067708
## 68 1239    68        0.07901282
## 69 1671    69        0.12769372
## 70  764    70        0.07724194
## 71 1311    71        0.13732062
## 72 1539    72        0.10431777
## 73 1297    73        0.07029810
## 74 1476    74        0.07464725
## 75 1853    75        0.10081062
## 76  724    76        0.04745674
## 77  925    77        0.08413680
## 78 1463    78        0.11211587
## 79 1241    79        0.07960742
## 80  881    80        0.04165879
## 81 1409    81        0.07333195
## 82  678    82        0.03242623
## 83  889    83        0.04047348
## 84  707    84        0.04426774
## 85  738    85        0.04636261
## 86  550    86        0.03138553
## 87 1026    87        0.02997196
## 88  789    88        0.02314054
## 89  654    89        0.02288474

write.csv(big_b_a_single_traj_data,
          "../Generated_Data/Representative_Simulations/big_b_a_single_sim_traj_data.csv",
          row.names = FALSE)

Load small b_a parameter combination

small_b_a_MLE = read.csv("../Generated_Data/Profiles/N_12_Model/top_2_LL_data/Man_Table_data/b_a_profile_top_2_LL_via_case_and_antibody_LL_small_b_a_parameter_combination_for_simulation.csv")
small_b_a_MLE

##   X        b_a M_0 V_0 K_0      R_0       b_q      b_p      p_S p_H_cond_S
## 1 1 0.06896552   5  13  14 6.098822 0.2343746 0.940249 0.148861  0.1569547
##   phi_E phi_U phi_S   h_V    gamma   N_0      E_0      z_0 C_0
## 1  1.09  1.09   0.2 0.125 6.333218 8e+06 63566.34 13443.03   0
##   social_distancing_start_time quarantine_start_time PCR_sens   sigma_M
## 1                           17                    22      0.9 0.2709327
##     beta_w_3  beta_w_2  beta_w_1 beta_w_0    g_0       g_F sigma_epsilon
## 1 0.01215073 0.9810086 -37.23481 229.4094 1183.3 0.1162005      109.1121
##   G_w_y_scaling  msg iter_num param_index  loglik nfail trace_num
## 1         0.162 mif1        2          44 -627.43    NA        NA
##    loglist.se Antibody_Mean_LL Antibody_LL_SE Median_Herd_Immunity
## 1 0.008839892        -24.41689    0.002390176            0.1998329
##   combo_num sim_subset_index duration_of_symp_1 duration_of_symp_2
## 1        13                1                  5          0.1578976
##   duration_of_symp gamma_total     Beta    R_0_P    R_0_A  R_0_S_1
## 1         5.157898   0.1938774 1.182424 1.019975 0.347037 0.880084
##      R_0_S_2  R_0_NGM
## 1 0.02343045 2.270527

small_b_a_MLE_comb =   small_b_a_MLE %>%
  dplyr::select(-msg, -iter_num, -param_index,
                -loglik, -nfail, -trace_num, -loglist.se,
                -Antibody_Mean_LL, -Antibody_LL_SE,
                -Median_Herd_Immunity, -combo_num,
                -sim_subset_index, -duration_of_symp_1,
                -duration_of_symp_2, -duration_of_symp,
                -gamma_total, -Beta,
                -R_0_P, -R_0_A, -R_0_S_1, -R_0_S_2,
                -R_0_NGM)

Simulate from small b_a parameter combination

##Simulation from ML
sim_data_small_b_a = simulate(nsim = 101,
                    seed = 12345,
                    times = Observed_data$times,
                    t0 = t0,
                    rprocess = pomp::euler(rproc,delta.t = 1),
                    params = small_b_a_MLE_comb,
                    paramnames = paramnames,
                    statenames = statenames,
                    obsnames = obsnames,
                    accumvars = acumvarnames,
                    rinit = init,
                    rmeas = rmeas,
                    partrans = par_trans,
                    covar = covar,
                    format = "data.frame")

sim_data_small_b_a_traj_data = sim_data_small_b_a %>%
  dplyr::select(time, .id, Y) 

p = ggplot(data = sim_data_small_b_a_traj_data,
           aes(x = time,
               y = Y,
               color = .id)) +
  geom_point() + geom_line(aes(group = .id)) +
  geom_point(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  geom_line(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  rahul_man_figure_theme +
  theme(legend.position = "None")
p

png("../Figures/Representative_Simulations/small_b_a_all_traj_sim.png")
print(p)
dev.off()

## quartz_off_screen 
##                 2

Plot 6 individual trajectories from small b_a parameter combination

select_trajectories  = filter(sim_data_small_b_a_traj_data, .id %in% seq(from = 5, to = 10))
select_trajectories = dplyr::select(select_trajectories, time, .id, Y)
select_trajectories$type = "Sim"

library(RColorBrewer)
full_blue_pallete = brewer.pal(9, "Blues")
sim_traj_pallete = full_blue_pallete[9:4]
sim_traj_scale = scale_color_manual(name = "Legend", values = c( sim_traj_pallete), labels =  c("Sim_Traj_1", "Sim_Traj_2", "Sim_Traj_3", "Sim_Traj_4", "Sim_Traj_5", "Sim_Traj_6")) 
p = ggplot(data = select_trajectories,
           aes(x = time, y = Y, color = .id)) + geom_point(size = 2) + geom_line(aes(group = .id)) + rahul_theme + rahul_man_figure_theme + theme_white_background +
  sim_traj_scale
p

png("../Figures/Representative_Simulations/small_b_a_5_traj_sim.png")
print(p)
dev.off()

## quartz_off_screen 
##                 2

p = ggplot(data = select_trajectories,
           aes(x = time, y = Y, color = .id)) + geom_point(size = 2) + geom_line(aes(group = .id)) + rahul_theme +
  rahul_man_figure_theme + theme_white_background +
  sim_traj_scale + geom_point(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  geom_line(data = Observed_data,
                              aes(x = times, y = Y), color = 'red')
p

png("../Figures/Representative_Simulations/small_b_a_5_traj_sim_with_obs.png")
print(p)
dev.off()

## quartz_off_screen 
##                 2

p = ggplot(data = select_trajectories,
           aes(x = time, y = Y, color = .id)) + geom_point(size = 2) + geom_line(aes(group = .id))  +
  theme_white_background +
  sim_traj_scale +
    facet_wrap(~.id, ncol = 1) +
  geom_point(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  geom_line(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  theme(legend.position = "None") +
  theme(strip.background = element_blank(),
        strip.text.x = element_blank())
p

png("../Figures/Representative_Simulations/small_b_a_5_traj_sim_with_obs_facet.png")
print(p)
dev.off()

## quartz_off_screen 
##                 2

p = ggplot(data = select_trajectories,
           aes(x = time, y = Y, color = .id)) + geom_point(size = 2) + geom_line(aes(group = .id))  +
  theme_white_background +
  sim_traj_scale +
    facet_wrap(~.id, ncol = 1) +
  geom_point(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  geom_line(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  theme(legend.position = "None") 
p

Pick “good” simluation

(One that corresponds well with the data by eye)

small_b_a_single_traj = select_trajectories %>%
  filter(.id == "5")
p = ggplot(data = small_b_a_single_traj,
           aes(x = time, y = Y)) + geom_point(size = 2,
                                              color = 'blue') +
  geom_line(color = 'blue') + rahul_man_figure_theme +
  geom_point(data = Observed_data,
                              aes(x = times, y = Y), color = 'red') +
  geom_line(data = Observed_data,
                              aes(x = times, y = Y), color = 'red')
p

png("../Figures/Representative_Simulations/small_b_a_best_traj_vs_obs.png")
print(p)
dev.off()

## quartz_off_screen 
##                 2

Save single trajectory for high b_a parameter combination

small_b_a_single_traj_data = small_b_a_single_traj %>%
  dplyr::select(Y =Y, times = time) %>%
  join(nyc_rel_testing_data) %>%
  mutate(obs_prop_positive = Y/Total_Number_of_Tests_Performed) %>%
  dplyr::select(-Total_Number_of_Tests_Performed)

## Joining by: times

small_b_a_single_traj_data$obs_prop_positive[1] = NA
small_b_a_single_traj_data

##       Y times obs_prop_positive
## 1     0     1                NA
## 2     0     2        0.00000000
## 3     3     3        0.37500000
## 4     6     4        0.16666667
## 5     1     5        0.02222222
## 6    12     6        0.18750000
## 7     4     7        0.06349206
## 8     7     8        0.05384615
## 9     7     9        0.03954802
## 10   34    10        0.19540230
## 11   28    11        0.07547170
## 12   94    12        0.13603473
## 13  173    13        0.33855186
## 14  179    14        0.21723301
## 15  329    15        0.36924804
## 16  476    16        0.21597096
## 17 1273    17        0.38068182
## 18 2050    18        0.42320396
## 19 2388    19        0.47035651
## 20 1435    20        0.20716039
## 21 1560    21        0.23586332
## 22 1732    22        0.32151476
## 23 3112    23        0.54730918
## 24 4744    24        0.68209921
## 25 2696    25        0.38942655
## 26 2704    26        0.34021137
## 27 5900    27        0.81921688
## 28 5433    28        0.87614901
## 29 3917    29        0.39597655
## 30 6875    30        0.93779839
## 31 6490    31        0.86383602
## 32 4562    32        0.48723700
## 33 5765    33        0.51253556
## 34 4623    34        0.59163041
## 35 5977    35        0.73490717
## 36 4325    36        0.48275477
## 37 7008    37        0.66288309
## 38 3601    38        0.33666791
## 39 5209    39        0.47238596
## 40 6349    40        0.53978915
## 41 2364    41        0.21923398
## 42 3607    42        0.44790761
## 43 4012    43        0.46061998
## 44 3257    44        0.20092535
## 45 6057    45        0.48363143
## 46 4975    46        0.42318816
## 47 3695    47        0.31714016
## 48 3289    48        0.32022198
## 49 2584    49        0.33136702
## 50 3370    50        0.38169668
## 51 3098    51        0.27525544
## 52 3735    52        0.27948219
## 53 2423    53        0.13375656
## 54 2157    54        0.12513053
## 55 2300    55        0.16956650
## 56 2525    56        0.23954084
## 57 1919    57        0.19399515
## 58 1848    58        0.13137129
## 59 2402    59        0.16841958
## 60 1605    60        0.11894175
## 61 1658    61        0.10077802
## 62 1585    62        0.11864661
## 63 1925    63        0.20962648
## 64 2196    64        0.20301377
## 65 2624    65        0.20113445
## 66 1331    66        0.08131224
## 67 2298    67        0.14960937
## 68 1285    68        0.08194630
## 69 1897    69        0.14496408
## 70 1399    70        0.14144171
## 71 1147    71        0.12014245
## 72 1507    72        0.10214872
## 73  757    73        0.04102981
## 74 1312    74        0.06635311
## 75 1262    75        0.06865785
## 76  910    76        0.05964866
## 77 1007    77        0.09159542
## 78 1260    78        0.09655912
## 79 1240    79        0.07954327
## 80  814    80        0.03849064
## 81  862    81        0.04486312
## 82  920    82        0.04400019
## 83 1038    83        0.04725700
## 84  662    84        0.04145013
## 85  719    85        0.04516899
## 86 1054    86        0.06014609
## 87  462    87        0.01349614
## 88  499    88        0.01463515
## 89  779    89        0.02725873

write.csv(small_b_a_single_traj_data,
          "../Generated_Data/Representative_Simulations/small_b_a_single_sim_traj_data.csv",
          row.names = FALSE)